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Ground-state phase diagram of the square lattice Hubbard model from density matrix embedding 

theory 
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We compute the ground-state phase diagram of the Hubbard and frustrated Hubbard models on the square 
lattice with density matrix embedding theory using clusters of up to 16 sites. We provide an error model to 
estimate the reliability of the computations and complexity of the physics at different points in the diagram. 
We find superconductivity in the ground-state as well as competition between inhomogeneous charge, spin, and 
pairing states at low doping. The estimated errors in the study are below Tc in the cuprates and on the scale of 
contributions in real materials that are neglected in the Hubbard model. 


The Hubbard model OHD is one of the simplest quan¬ 
tum lattice models of correlated electron materials. Its one- 
band realization on the square lattice plays a central role 
in understanding the essential physics of high temperature 
superconductivity (H 13 . Rigorous, near exact results are 
available in certain limits (61: at high temperatures from 
series expansions CHm, in infinite dimensions from con¬ 
verged dynamical mean-field theory (TJHMI, and at weak 
coupling from perturbation theory ca and renormalization 
group analysis ||l6l[I3. Further, at half-filling, the model has 
no fermion sign problem, and unbiased determinantal quan¬ 
tum Monte Carlo simulations can be converged M- Away 
from these limits, however, approximations are necessary. 
Many numerical methods have been applied to the model at 
both finite and zero temperature, including fixed node, con¬ 
strained path, determinantal, and variational quantum Monte 
Carlo (QMC) (1914291 , density matrix renormalization group 
(DMRG) (231321, and dynamical cluster (DCA) (33l[3l and 
(cluster) dynamical mean-field theories (CDMFT) (35l l36l . 
These have revealed rich phenomenology in the phase dia¬ 
gram including metallic, antiferromagnetic, d-wave (and other 
kinds of) superconducting phases, a pseudogap regime, and 
inhomogeneous orders such as stripes, and charge, spin, and 
pair-density waves ID 

Here, we employ density matrix embedding theory 
(DMET) (J7l l38l . together with clusters of up to 16 sites 
and thermodynamic extrapolation, to compute a calibrated 
ground-state phase diagram for the Hubbard model. We use 
the term calibrated as we provide an error model to estimate 
the quality of the results, and by proxy, the complexity of the 
underlying physics. The one-band (frustrated) Hubbard model 
on the L X L square lattice is 

{ij)a {{ij))cr i 

where (...) and ((...)) denote nearest and next-nearest neigh¬ 
bors, respectively, a- ^ destroys (creates) a particle on site i 
with spin cr, and n^cr = is the number operator. The 

standard Hubbard model corresponds to f' = 0 (we fix t = 1). 
We further study the frustrated model with t' = ±0.2. 

DMET is a cluster impurity method which is exact for 
weak coupling (U = 0) and weak hybridization (t = 0) 
and which becomes exact for arbitrary U as the cluster size 


Nc increases. It differs from Green function impurity meth¬ 
ods such as the DCA or (C)DMFT, as it is a wavefunction 
method, with a finite bath constructed to reproduce the entan¬ 
glement of the cluster with the remaining lattice sites without 
bath discretization error. DMET has recently been applied 
and benchmarked in a variety of settings from lattice mod¬ 
els OUllSID to ab-initio calculations with realistic long- 
range interactions isaiia, and for ground-state and spectral 
quantities m. In its ground-state formulation, the use of 
wavefunctions substantially lowers the cost relative to Green 
function impurity methods, allowing larger clusters to become 
computationally affordable. 

We briefly summarize the method here, with details in the 
supplementary information and original references (J7l |38]| . 
DMET maps the problem of solving for the bulk ground-state 
|T^) (on the L X L lattice for L sufficiently large) to solving for 
the ground-state of an impurity model with Nc impurity and 
Nc bath sites. The exact mapping is defined via the Schmidt 
decomposition HU of the exact |T^) = \ \^i) \^i)^ where 

{\ai)} denotes impurity states, and {\bi)}, bath states. The 
bulk |T^) can be expressed exactly in the Schmidt subspace 
{\ciibj)} and is the ground state of the impurity Hamiltonian 
defined as i^imp = PHP, P = thus es¬ 

tablishing the exact ground-state bulk to impurity mapping. 
In practice, however, the exact |T^) is, of course, unknown! 
DMET therefore solves an approximate impurity problem de¬ 
fined from a model bulk wavefunction |4>), the ground-state 
of a quadratic Hamiltonian h = where Hq is one-body 

part of the Hubbard Hamiltonian, and is a one-body operator 
acting in each cluster unit cell of the bulk lattice, to be deter¬ 
mined. Via |4>) we define the bath space, impurity Hamilto¬ 
nian, and impurity model ground-state |T^') (which is now an 
approximation to the exact bulk wavefunction |T^)) and from 
which energies and local observables can be measured. Under 
this approximation, the bath Hilbert space spanned by 
(of equal size to the impurity Hilbert space) becomes isomor¬ 
phic to the Eock space of Nc (one-particle) sites, i.e. the bath 
sites. All these quantities are functions of the one-body op¬ 
erator u, which is determined self-consistently by matching 
the one-particle density matrix of the impurity wavefunction 
|T^'(i4)), and the model lattice wavefunction |4>(i4)), corre¬ 
sponding to the optimization min^^ \ (i^) |at Oj |T^(i4)) — 

{^{u)\alaj\^{u)) p, where i, j label impurity and bath sites. 
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In this work, we used two small modifications of the original 
DMET procedure in Ref. IIJTII . First, we allowed u to vary 
over particle non-conserving terms, thus allowing |^(r^)) to 
spontaneously break particle number symmetry in order to de¬ 
scribe superconducting phases. Second, we used an additional 
chemical potential on the impurity sites, to ensure that the im¬ 
purity fillings for 1^) and |^') exactly match. 

To obtain the ground-state phase diagram, we carried out 
DMET calculations using 2x2, 4x2, 8x2, and 4x4 impurity 
clusters, cut from a bulk square lattice with L = 72. We con¬ 
sidered t' = 0, 0.2, —0.2, and U = 2,4, 6, 8, and various den¬ 
sities between n = 0.6 — 1. The impurity model ground-state 
|T^') was determined using a DMRG solver with a maximum 
number of renormalized states M = 2000, and which allowed 
for U{1) and SU{2) symmetry breaking. The energy, local 
moment m = — n^), double occupancy D = 

and local d-wave pairing dsc = 
measured from |T^'). 

The finite cluster DMET energies and measurements con¬ 
tain 3 sources of error relative to the exact thermodynamic 
limit. These are from (i) errors in the DMET self-consistency, 
(ii) finite M in the DMRG solver (only significant for the 
8x2 and 4x4 clusters, corresponding to 32 impurity plus 
bath sites in the impurity model), which also induces error 
in the correlation potential u, (iii) finite impurity cluster size. 
(The error from the use of a finite 72 x 72 bulk lattice, is so 
small as to not affect any of the significant digits presented 
here). To estimate the thermodynamic result, we (i) estimated 
DMET self-consistency quality by the convergence of expec¬ 
tation values in the last iterations, (ii) extrapolated DMRG en¬ 
ergies and expectation values at finite M to infinite M, using 
the linear relation with DMRG density matrix truncation er¬ 
ror 1461 , (iii) estimated the error in u due to finite M, by ex¬ 
trapolating expectation values from self-consistent u{M) ob¬ 
tained with different solver accuracy, (iv) extrapolated cluster 
size to infinite size, with the scaling Nc appropriate to a 
non-translationally-invariant impurity. Each of (i) to (iv) gives 
an estimate of an uncertainty component (for linear extrapola¬ 
tions, we use the la standard deviation), which we combined 
to obtain a single error bar on the DMET thermodynamic es¬ 
timates. Details of the error estimation and a discussion of the 
complete data (of which only a fraction is presented here) are 
given in the supplementary information. 

We first verify the accuracy of our thermodynamic esti¬ 
mates and error bars by comparing to benchmark data avail¬ 
able at half-filling. In Table and Fig. we compare the 
DMET ground-state energy, double occupancies, and stag¬ 
gered magnetization with exact estimates at half-filling, as ob¬ 
tained from ground-state (auxiliary field) determinantal QMC 
(AFQMC) calculations on finite square lattices extrapolated 
to infinite size (13, and DMRG on long open cylinders, ex¬ 
trapolated to infinite width and length 1491 . For comparison, 
we also show recent DCA energies computed at the lowest 
published temperatures, T = 0.05 — 0.15t l5Ql . 

The data shows the high accuracy of the DMET energies at 


TABLE I. Ground state energy of the 2D Hubbard model. All the 
numbers are extrapolated to the thermodynamic limit. (CP-)AFQMC 
results are from Zhang w\ . Note that the half-filling results do not 
involve the constrained path approximation l48l , thus is numerically 
exact. DMRG results are from White 1491 . 


U/t 

Filling 

DMET AFQMC 

CP-AFQMC 

DMRG 

2 

1.0 

-1.1764(3) -1.1763(2) 

- 

-1.176(2) 

4 

1.0 

-0.8604(3) -0.8603(2) 

- 

-0.862(2) 

6 

1.0 

-0.6561(5) -0.6568(3) 

- 

-0.658(1) 

8 

1.0 

-0.5234(10) -0.5247(2) 

- 

-0.5248(2) 

12 

1.0 

-0.3686(10) -0.3693(2) 

- 

-0.3696(3) 

4 

0.8 

-1.108(2) 

-1.110(3) 

-1.1040(14) 

4 

0.6 

-1.1846(5) 

-1.185(1) 

- 

4 

0.3 

-0.8800(3) 

-0.879(1) 

- 



(b) (c) 


FIG. 1. (color online) Benchmark and uncertainties for t' — Hub¬ 
bard model, (a) Energy at half-filling. Ground state estimates from 
DMET, AFQMC l47l and DMRG |4^ , compared to a recent DCA 
study Bol . The temperatures are the lowest published values in the 
DCA study. (DCA data at U=8 is from 50-site impurity cluster cal¬ 
culations, not extrapolated to thermodynamic limit.) (b) Energy un¬ 
certainties across the phase space. The areas of the circles are pro¬ 
portional to the estimated uncertainties, (c) Staggered magnetization 
(m) and double occupancy(i9) at half-filling. The solid blue line is 
the spin-1 Heisenberg limit m = 0.3070(3) 1^ . 

half-filling. The error bars from DMET, AFQMC, and DMRG 
are all consistent. Except for 17 = 8 where the error is slightly 
larger, DMET provides the same number of significant dig¬ 
its as the “exact” AFQMC number with an accuracy better 
than O.OOIL As a point of reference, the uncertainties in the 
ground-state methods are significantly smaller than the finite 
temperature contributions to the low-temperature DCA calcu¬ 
lations (Fig. [Ja)). 

Figure ^c) further gives the half-filling staggered magne¬ 
tization and double occupancies computed with DMET, as 
compared with AFQMC. The DMET double occupancies are 
obtained with similar error bars to the “exact” AFQMC es¬ 
timates. The DMET staggered magnetization, a non-local 
quantity, exhibits larger errors at the smallest 17 = 2 (a cluster 
size effect) but for 17 >4 appears similarly, or in fact more 
accurate than the AFQMC result. At the largest value 17 = 12, 
we find m = 0.327(15), slightly above the exact Heisenberg 
value (53 

The half-filling benchmarks lend confidence to the DMET 
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thermodynamic estimates of the energy and observables, and 
their associated error bars. We therefore use the same error 
model to estimate the accuracy of the DMET energies and 
expectation values away from half-filling, in the absence of 
benchmark data. Although exact thermodynamic limit re¬ 
sults are not available away from half-filling, we can verify 
our error model at low density using constrained path (CP) 
AFQMC, a sign-free QMC with a bias that disappears at low 
density and small U |[23l|24l. For U = A and n < 0.6, a 
parameter regime where CP-AFQMC is very accurate, the 
DMET and CP-AFQMC energies agree to O.OOlt (Table |I]). 
Fig. [^b) shows the energy uncertainties across the phase di¬ 
agram for t' = 0. (The same figure for t' = ±0.2 is given 
in the supplementary information, which, in general, displays 
smaller error than t' = 0). As expected, the accuracy away 
from half-filling is significantly lower than at half-filling, with 
the largest errors found in the underdoped region of n=0.8-0.9. 
The main source of error is from cluster size extrapolation, es¬ 
pecially in the underdoped region. Farge errors can be viewed 
as reflecting underlying physics, as they coincide (see below) 
with phase boundaries and/or the onset of competing inhomo¬ 
geneous orders, both of which are sensitive to cluster shape, 
and thus lead to errors in extrapolation. 
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FIG. 2. (color online) Phase diagrams of the standard and frus¬ 
trated Hubbard model. Orders are represented with three primary 
colors: red (antiferromagnetism), green (d-wave superconductivity) 
and blue (inhomogeneity), with the brightness proportional to the ro¬ 
bustness of the order (discussed in the supplementary information). 
The points highlighted with letters: (a) local phase separation; (b) 
d-wave SC with a slight modulation in (tt, tt) direction; (c) SC with 
a weak spin density wave (SDW); (d) a “classic” stripe phase; (e) 
stripe with pair-density wave (PDW) coexisting with SC; (f) CDW 
and spin 7r-phase shift; (g) intermediate points between AF and SC 
where both order parameters extrapolate to zero. 

We present the DMET phase diagrams in Fig. We ob¬ 
serve (i) an AF phase at half-filling, (ii) a metallic phase at 
large dopings and at small U, enhanced by frustration, (iii) a 
region of d-wave SC order at intermediate dopings and suf¬ 
ficiently large U, (iv) a region of coexisting AF and SC or¬ 
der, (v) a region with various inhomogeneous charge, spin, 
and superconducting orders, (vi) points in between the AF 
and SC phase where the AF and SC orders extrapolate to 
zero. (The metallic phase is predicted, to be unstable at very 
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FIG. 3. (color online) Antiferromagnetic (red circle) and (d-wave) 
superconducting (green square) order parameters at U=4. 


weak coupling and large dopings from weak coupling expan¬ 
sions Oil El, but this is associated with an exponentially 
small energy scale not probed here). Fig. shows the average 
AF and d-wave SC order parameters as a function of filling for 
U=A. We find that t' = 0.2 stabilizes AF versus SC, and the 
reverse is true for t' = —0.2. For t' = 0, the peak in SC order 
is around (n) = 0.9 and SC extends to (n) ^0.8. The figures 
also show the suppression (enhancement) of SC order with 
positive (negative) t'. As positive (negative) t' corresponds 
to electron-(hole-)doped cuprates, our results are consistent 
with the stronger superconductivity found in hole-doped ma¬ 
terials |[53ti55l . 

The presence of SC in the Hubbard model ground-state has 
previously been much discussed. From the Mermin-Wagner 
theorem, long-range order is not allowed at finite tempera¬ 
tures |[56H58]| . but at zero temperature, such long-range or¬ 
der can exist. In a cluster mean-field approach embodied by 
cluster DMET (and similarly CDMFT) a concern is that the 
observation of local order in finite clusters does not translate 
into true long-range order. However, our estimates indicate 
that a homogeneous SC order parameter survives in the infi¬ 
nite cluster limit, within the error bars of our extrapolation. 

At t' = 0, we observe a banana-shaped SC region. At 
U = 6 and n = 0.875 (between the AF and SC phases) we 
find that the AF and SC order parameters are nonzero in fi¬ 
nite clusters, but extrapolate to 0 in the thermodynamic limit. 
However, for the analogous U = 8, n = 0.875 point, a SC 
state with strong inhomogeneity appears which creates large 
uncertainties in the extrapolated order parameters, thus the 
precise location of the SC phase boundary at f/ = 8 is un¬ 
certain. 

We now further discuss the intermediate region between the 
AF and SC phases (low doping and large U). In this region, 
a variety of spin-density (251 l59l463]| charge-density (251164]- 
l66ll . pair-density wave (66ti69ll , and stripe orders (^[^[63l 
|701473]| , have been posited in both the Hubbard model and the 
simpler t-J model. These inhomogeneous phases are pro¬ 
posed to be relevant in the pseudogap physics i67i Ea m- 
l78l . Recent projected entangled pair state (PEPS) studies of 
the t-J model suggest that different inhomogeneous and ho¬ 
mogeneous states are near degenerate at low doping ca. Our 
work indicates that the Hubbard model behaves similarly. For 
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FIG. 4. (color online) Local order parameters in the (frustrated) Hub¬ 
bard model at selected points in the strong coupling regime (U =8). 
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large U and low doping n = 0.875 — 0.8we find some points 
(marked (g) in Fig.[^ where the AF and SC order parameters 
are small or vanish, but also many other points with various 
inhomogeneous orders. Some representative examples of in¬ 
homogeneous states are shown in Fig. These correspond 
to (i) a local phase separation between a half-filled, antiferro¬ 
magnetic phase and a superconducting ribbon (Fig.Qa)), (ii) a 
classic stripe phase order, with a period 4 charge and period 8 
spin density modulation (Fig.Hlb)), very similar to as seen in 
earlier DMRG ladder studies . There is also a coexisting 
weak PDW (exhibiting a sign change across the cell), consis¬ 
tent with earlier stripe proposals (691. (iii) Inhomogeneities 
in the pairing order coexisting with the charge and spin or¬ 
ders in, eg., Fig.|^c) where we see a PDW with an 8 unit cell 
wavelength coexisting with a CDW with a 4 unit cell wave¬ 
length and a 8 unit cell SDW. The PDW is all positive (on the 
ladders) indicating coexistence with superconductivity, simi¬ 
lar to a recent theoretical proposal (see e.g. Ref. (681). The 
inhomogeneity is mainly observed with zero or negative t', 
corresponding to the hole-doped cuprates. Fig. [Jd) shows 
an example with t' = 0.2, where the inhomogeneity is much 
weaker than in the t' < 0 cases. Although only 8x2 clusters 
are shown above, not all 8 x 2 clusters are inhomogeneous, and 
similarly not all 4 x 4 clusters are homogeneous. A detailed 
analysis of observed inhomogeneous phases, and the determi¬ 
nation of the phase diagram, is presented in the supplementary 
information. While the impurity clusters we use are still too 
small to definitively resolve the competing orders, they hint at 
the possible behaviour and energy resolution required to de¬ 
termine the ground state at various points in the phase space, 
and where we should focus our attention using larger clusters 


in future studies. 

To summarize, we have computed a ground-state phase di¬ 
agram for the Hubbard and frustrated Hubbard models on the 
square lattice using cluster DMET. At half-filling, the accu¬ 
racy achieved by DMET appears competitive with the best ex¬ 
act benchmarks, while away from half-filling our error model 
suggests that the calculations remain very accurate. We ob¬ 
serve standard AE and metallic phases, regions of d-wave SC 
pairing order, and several kinds of inhomogeneities. At spe¬ 
cial points in the phase diagram, the inhomogeneous and ho¬ 
mogeneous solutions associated with 8x2 and 4x4 clusters 
are very close in energy and definitive characterization will re¬ 
quire higher energy resolution with larger clusters. However, 
for real materials such as the cuprates, assuming t ^ 3000K, 
the energy resolution achieved here for most of the phase dia¬ 
gram is already below the superconducting transition temper¬ 
ature, suggesting that the near degeneracy of these orders will 
be lifted by terms beyond those in the Hubbard model, such 
as long range charge and hopping terms, multi-orbital effects, 
and interlayer coupling. Moving beyond the Hubbard model 
to more realistic material models thus now appears of princi¬ 
pal relevance. 
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I. SUMMARY OF DMET 

Fig. 1 illustrates the computational flow of a DMET calcu¬ 
lation. A DMET self-consistency cycle consists of (i) solv¬ 
ing for the ground-state of the DMET lattice Hamiltonian, (ii) 
building the impurity Hamiltonian, (iii) solving for the impu¬ 
rity Hamiltonian ground-state and observables, and (iv) fitting 
the DMET correlation potential. As discussed in the main 
text, in this work we allow the DMET solutions to sponta¬ 
neously break particle number and spin symmetry. We also 
include a chemical potential in the self-consistency. Here, we 
explain some general aspects of practical DMET calculations 
which have not been discussed in detail in the existing liter¬ 
ature, as well as describe the technical extensions to broken 
particle number symmetry, and the self-consistency procedure 
for the additional chemical potential. 



FIG. 1. Graphical representation of the DMET procedure. 


A. DMET correlation potential 

A general DMET correlation potential is a quadratic oper¬ 
ator. It is local in the sense that it does not have cross terms 
between different images of the impurity cluster on the lattice. 
In the original DMET paper^, it took the form 

U = Y,UC = H (1) 

C C i,jec,a 

In Eq. (I), C ranges over all impurity cluster supercells within 
the (large) lattice, i,j range over sites in the same cluster C, 
and a G denotes the two flavors of spin. In this form, 

the correlation potential has Nc{Nc + l)/2 free parameters 


(here and later on, we assume real potentials) where Nc is the 
number of impurity cluster sites. Eor spontaneously broken 
particle number and spin symmetry, the correlation potential 
acquires additional terms, 

U = + /l.C. (2) 

In this work, we only allow singlet pairing (strictly speaking, 
Sz = 0 pairing) but it is straightforward to extend the above 
to triplet pairing. The normal part v has two spin components. 
The pairing term A has free parameters (it is symmet¬ 
ric when spin symmetry is preserved, but we allow for spin 
symmetry breaking). In total, the correlation potential u has 
Nc{2Nc + 1) degrees of freedom. 


B. DMET lattice Hamiltonian 


The DMET lattice Hamiltonian (including a chemical po¬ 
tential /in) is 

h h U /in ^ ^ (3 tz.c. (3) 

ija 


where h is the hopping term in the (frustrated) Hubbard 
model, h' can be rewritten in the form of of a spin-unrestricted 
Bogoliubov-de Gennes (BdG)^’^ equation. 



These coupled equations are expressed concisely as 




Ua K. 
V^/3 


Ua V^\ fSa \ 

Vp J\ -^/3 J 

(5) 


where 5c, and are both positive, h' is diagonalized by trans¬ 
forming to the Bogoliubov quasiparticles. 




= u3a]p + 


( 6 ) 


Note that the number of {cj^} and {c^} quasiparticles will 
differ if 5'^ 7 ^ 0 in the physical ground-state. 

In terms of the quasiparticles, the lattice Hamiltonian in 
Eq. (3) is diagonalized as 


h — Eq H- ^ ^ ^i<j^i(j^i(T 
i(j 


( 7 ) 
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and the (ground state) quasiparticle vacuum |—), defined by 
Cicr\—) = 0, has energy Eq. The quasiparticle vacuum is 
also known as the Bardeen-Cooper-Schrieffer (BCS) ground- 
state"^. 


C. DMET impurity model Schmidt subspace 


A and B can also be obtained directly from the one-particle 
density matrix. The rotation between C and Co leaves the 
one-body density matrix invariant, thus 

p = (ajaj) = = CC^ 

^(AA^ AB^ A /^An,p pM (12) 

( BA^ BB'^ + DD^ ) ~\ Pc Penv ) 


We now discuss how to define the impurity model Schmidt 
subspace corresponding to a BCS ground-state of the lat¬ 
tice Hamiltonian in Eq. (3). To start, we review the “prod¬ 
uct space” construction of the impurity model Schmidt sub¬ 
space, starting from the lattice Hamiltonian Slater determinant 
ground-state, as used in the original DMET^’^. 

The original DMET impurity model consists of a set of im¬ 
purity sites augmented by a set of bath modes. In Ref.^, the 
bath modes are defined through the projected overlap matrix 
of the Slater determinant. We compute the projected overlap 
matrix from the Slater determinant coefficient matrix, 

c.= (") (8) 

where the rows denote physical sites, and columns are occu¬ 
pied modes (orbitals). The upper part M has Nc rows. The 
projected overlap matrix is 


S = 


(9) 


Erom the singular value decomposition (SVD) of M as M = 
LER^ (where we use the “full” form of the SVD, L is A^c x 
Nc, E is Nc X n, and Ris n x n) then S = R{'E^'E)R^, i.e 
R is the eigenvector matrix of the projected overlap matrix. R 
defines a unitary transformation of the occupied modes in Co, 
giving a new coefficient matrix C = CqR, where 


Defining the eigendecomposition pimp = UAU^, we find 

A = UA^ and B = Pc{A^)~^ (13) 


The above defines the impurity model Schmidt subspace as 
a tensor product of the impurity site space, and a bath space, 
thus we refer to it as a “product-space” embedding construc¬ 
tion. However, for the BCS state, it is easier to use a slightly 
different, but equivalent construction. We explain this first for 
the Slater determinant. Here we build a L x 2Nc matrix C 2 , 
whose columns span the same vector space as Ci in Eq. (11), 
but which does not have the block structure. We start with the 
‘hole” one-particle density matrix 

Ph = [aitt]) = I - = I - p (14) 

We can replace p with ph in Eqs. (12) and (13) and compute 
an analogous set of coefficients A' and B'. Taking A, B, and 
A', B' gives C 2 , 


C 2 


A A' \ 


(15) 


The 2Nc columns of C 2 span exactly the same space as Ci 
(proved in the Appendix). Thus, we can equivalently define 
the Schmidt subspace from the columns of C 2 as we can from 
Cl. Transforming to the quasiparticle vacuum of the Slater 
determinant, | —), the columns of C 2 define a set of 2Nc quasi¬ 
particle creation operators 


C = 


( 1 

A 

0 \ 

1 NR 

-[b 

D J 


( 10 ) 


and the second equality follows because E is a rectangu¬ 
lar matrix of the form (diag(cr), 0, 0,... 0), where the first 
Nc columns constitute a diagonal matrix, and the remaining 
n — Nc columns are zero columns. The first Nc columns 


of C, 


A 

B 


define the embedding modes, which have non¬ 


zero weight on the impurity sites. The matrix B defines the 
bath modes, which may be orthonormalized using the QR de¬ 
composition, B = QR. The remaining columns in C define 
the core modes, which have no weight on the impurity. The 
Schmidt subspace is then E{a\) ^ jF{h\) 0 |ei... Cu-nQ^ 
where {at} create the impurity site basis, {b]} create the bath 
modes (from the columns of Q), and |ei... Cu-nQ is the core 
state, defined by the columns of D. The coefficients defining 
{at}, {^t} can be gathered in the columns of a matrix Ci, 





^ji^ja' + Bjiajfj' 

(16) 

jEimp j£env 



(17) 

jEimp j£env 



that yields the Schmidt subspace as ^1“)- As 

the impurity model Schmidt subspace here does not (transpar¬ 
ently) separate between the impurity sites and environment 
sites, but rather involves a set of modes which are a linear 
transformation of both the occupied and virtual modes in the 
Slater determinant, we refer to this as a “quasiparticle em¬ 
bedding” construction. This provides an alternative view of 
the DMET embedding as an active space method that uses the 
embedding quasiparticles defined from C 2 as the active space, 
while freezing other excitations that involve only the environ¬ 
ment. 

Extending the quasiparticle embedding construction to 
BCS states is straightforward. By analogy with the one- 
particle density matrix of a Slater determinant, we define the 
generalized one-body density matrix for BCS states. 


CcT — 




1 Pa 


pa 


where is an Nc x Nc identity matrix. 


(18) 
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where the normal one-particle density matrices pa- = 
(al^aja) = VaVj = 1 — UaUj, and the pairing density 
matrix k = (ai^ajp) = i^a = = UaVj. The diagonal 

of G is formed by the hole and particle density matrices, and 
the off-diagonals by the pairing matrix. When the BCS state 
is a Slater determinant, /^ = 0. 

We reorganize the generalized density matrix G into impu¬ 
rity and environment blocks, placing the impurity (environ¬ 
ment) submatrices of p and n together. 

G= f (19) 

y Cjtc ^env / 


Then, similar to the treatment in Eq. (13), we rewrite the 
impurity part of the density matrix Gimp = AA^, and define 
a quasiparticle G 2 matrix 


C 2 


( Gcr,imp \ 
Yp',imp 

\ ^^,env 




A 

B 


2Lx2Nc 


( 20 ) 

where B = Gc{A^)~^. Eq. (20) defines a new set of 
quasiparticles (with associated quasiparticle creation opera¬ 
tors {c\^}) in Eq. (6) through the coefficients Ucj^Vcj. These 
are a unitary rotation of the original 2L quasiparticles such 
that only 2Nc of them have non-zero overlap with the impu¬ 
rity. As the rotation does not mix the quasiparticle creation 
and annihilation operators, the vacuum of Ci^j is still the BCS 
ground state |—). In analogy to the embedding for Slater de¬ 
terminants, the Schmidt subspace is now spanned by the em¬ 
bedding quasiparticles, 0 |—). 

To connect with Eq. (15), note that when the BCS state is a 
Slater determinant. Gimp and Gc are both block diagonal, and 
thus, A = dmg{A'^, A^),B = dmg{B'^,B^), and Eq. (20) 
becomes 


Ua = 

v-,= 


A'^ 0\ 
B'^ 0 ) 

0 A^ \ 
0 B^ J 


( 21 ) 

( 22 ) 


Combining both sets of spins, the quasiparticles in Eq. (22) 
then span exactly the same Hilbert space as the basis defined 
in Eq. (15). Eor general BCS ground states, however, A and B 
are not block diagonal, and the embedding quasiparticles are 
mixtures of particles and holes. 


D. DMET impurity Hamiltonian and DMRG solver 

Once the Schmidt subspace has been defined, the DMET 
Hamiltonian is formally obtained by projecting an interacting 
lattice Hamiltonian into the subspace as i^imp = PH'P, with 
the many-particle projector defined as 

P = ^\^nJ{^nJ, (23) 

P'icT 


where is a vector of occupation numbers of the embedding 
quasiparticles, and = Iln !“)• earlier 

DMET work, two choices of lattice Hamiltonian were used in 
the projection: the original interacting lattice Hamiltonian H 
(in this case the original Hubbard Hamiltonian) and a modified 
interacting lattice Hamiltonian H', where the interaction term 
U is only used in the impurity sites. As in earlier DMET work 
on lattice models, here we use the latter simpler Anderson- 
like lattice Hamiltonian H'. In H', on the environment sites 
(outside of the impurity cluster) the Coulomb interaction U is 
replaced by the correlation potential u, giving 


H' = h+ + '^ UriiaUip - /in (24) 

CT^imp iGimp 

The projection defined in Eq. (23) reduces to transforming 
{a-P i to the embedding quasiparticles using the inverse Bo- 
goliubov transformation, 

Act = uljcA + VijCj^ (25) 

and replacing the pure environment quasiparticle operators 
with their expectation values with the BCS ground state |—). 

After projection, we can write i^imp as a sum of one- and 
two-particle parts, iifimp = ^imp + ^imp, where /limp is 

^imp = T^ij^icj'-jcr + + ^0 (26) 


and c-p here denote the embedding quasiparticles. The two- 
particle part V contains many contributions due to the break¬ 
ing of particle number symmetry in the quasiparticle formula¬ 
tion. These have the form 


i^mp—2 ^ ^ '^pgsr,cr/xcJcr^g/x^s/xGcr 

pqsr,aiJ, 

P ^ ^ 3 “ Pi 




^ A A 

'^pqsr^pa^qa^spAp 


pqsr 




^pqsr,(J^pCF ^scr Gcr 




pq 


(27) 


Vimp connects N electron states with N^N ±2,N ±4 states. 
Eor brevity, we do not give the formula for the coefficients ex¬ 
plicitly (which are obtained by simple algebra from Eq. (25)). 
The scalar and one-particle terms in Vimp contain contributions 
from pure environment quasiparticles, and can be absorbed 
into /limp* 

We have adapted our quantum chemistry density matrix 
renormalization group (DMRG) code BLOCK^“^ to break 
U{1) particle number symmetry and to incorporate the Hamil¬ 
tonian terms in Eq. (26) and (27). While the full wavefunction 
is not restricted to U{1) symmetry, the particle quantum num¬ 
ber is still used in the calculation in the sense that the renor¬ 
malized states are required to carry a definite particle number. 
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This allows us to use the block-sparsity of the operators to 
tackle larger numbers of renormalized states. 

The quasi-particle basis associated with is not localized 
to a site, thus we use a localization and ordering procedure 
as used in quantum chemistry DMRG calculations to reduce 
long-range entanglement between the embedding quasiparti¬ 
cles. We find that the localization and reordering significantly 
reduce the DMRG truncation error, by up to a factor of 10. 


E. Expectation values 

As discussed in the original papers on DMET^, the DMET 
energy of i^imp defined in Eq. (26) does not correspond to the 
ground-state energy of the impurity cluster. This is because 
the impurity Hamiltonian contains three types of energy con¬ 
tributions: pure impurity, impurity-bath interactions, and pure 
bath (environment) parts. The proper DMET energy should 
exclude the pure environment contributions and include only 
part of the impurity-bath interaction energy. Therefore, the 
DMET energy is evaluated as a partial trace of the one- and 
two-particle reduced density matrices of the impurity wave- 
function. This partial trace can be equivalently implemented 
as a full trace, with appropriate scaling factors for terms in 
the Hamiltonian which couple the impurity and environment. 
Eor each class of term in the Hamiltonian, this scaling factor 
is given by the number of indices in the impurity, divided by 
the total number of indices. (Eor example, for the one-particle 
terms in the Hamiltonian, the contribution of the impurity- 
bath block to the total trace is scaled by a factor of |). 

An equivalent formulation for the Hubbard Hamiltonian 
(which contains no long-range Coulomb terms) is to evaluate 
the two-particle part of DMET energy as 

E2 = (^l^mpl^) = ^DMRG - (28) 

where |T^) is the DMRG ground state. Since /limp is a 
quadratic operator, E 2 can be computed only with knowledge 
of the DMRG energy and the one-particle (and pairing) den¬ 
sity matrix, avoiding explicitly evaluating (T^jf^mpl^) through 
the two-particle density matrix. 

The local spin moments and pairing are both one-particle 
quantities. We therefore obtain them from the one-particle 
and pairing density matrix p = {c\^Cja), of the 

DMRG wavefunction |T^), transformed back to the lattice site 
basis using Eq. (25). Note that p and n are defined not 

only for quasiparticles inside the impurity Schmidt subspace, 
but also for core quasiparticles. (In the quasiparticle approach, 
although p and n are themselves zero in the core, terms such 
as Qcj can appear in the expansion using Eq. (25) and result 
in non-zero expectation values due to the contributions of the 
new vacuum). If one is interested only in impurity cluster ex¬ 
pectation values, or for DMET lattice Hamiltonians without 
broken symmetry, the contribution of the core quasiparticles 
is strictly zero and may thus be omitted. However, for ordered 
(e.g. magnetic or superconducting) states, the core contribu¬ 
tion does not vanish an therefore cannot be neglected. Doing 
so would produce for example, the strange result of vanishing 


long-range correlations even in an long-range ordered DMET 
state. 

In this study, when a single value of the order parameter is 
given, it is computed using the 2 x 2 plaquette at the center 
of the impurity cluster, to minimize the boundary effects. The 
antiferromagnetic order parameter is defined as 

m = ^(mo,o + ^1,1 - ^0,1 - ^1,0) ( 29 ) 

and the d-wave parameter as 

d — ^[<^( 0 , 0 ),( 0 , 1 ) + <^( 1 , 0 ),( 1 , 1 ) ~ <^( 0 , 0 ),( 1 , 0 ) ~ <^( 0 , 1 ),( 1 , 1 )] 

(30) 

where rrii = 2(^.^ _ ^.^) and dij = -^{{oio^ajp) -h 
{o^jaCiip)) as defined in the main text. Of course in some cases 
there are also inhomogeneous states. When the inhomogene¬ 
ity is strong, we report here the full distribution of local order 
parameters. 

F. DMET self-consistency 

The DMET embedding constructs the impurity model via 
the model ground-state of the DMET lattice Hamiltonian, 
however, this state (and the lattice Hamiltonian) are functions 
of the correlation potential u. u is determined by the self- 
consistency procedure, which aims to minimize the difference 
between the embedding wavefunction and the DMET mean- 
field wavefunction, as measured by their (generalized) one- 
particle density matrix difference. In the quasiparticle embed¬ 
ding space, the one-particle and pairing density matrices of 
the mean-field wavefunction |T>) are simply zero. Concep¬ 
tually, the simplest technique is to define u to minimize the 
Erobenius norm, 

min||G'5,(„) - G$(„)||f 

= + 2|Kvi,,y p) 

However, as the derivative of the correlated wavefunction |T^) 
with respect to u is expensive, the above is solved in a two- 
step procedure consisting of an inner and outer loop. In the 
inner loop, we carry out min^^ ||G^ — G^(^u) IIf, i-e. the cor¬ 
related wavefunction is held fixed, while in the outer loop, the 
updated u leads to a new impurity model, and a new correlated 
wavefunction 

If the total particle number n is allowed to fiuctuate, as in a 
superconducting state, then one of the conjugate pairs (chem¬ 
ical potential) p or (particle density) (n) must be fixed. We 
usually want to express the observables as a function of dop¬ 
ing, or occupation, thus we fix (n) and determine the appro¬ 
priate p. Since the diagonal elements of the correlation and 
chemical potential appear redundant, how can one determine 
the chemical potential? Eormally, at the DMET mean-field 
level (Eq. (3)), there is a gauge freedom between u and /i, 
namely 

p'= 11 +(p,u'= u + (32) 
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however, this gauge freedom is lost at the embedding stage 
(Eq. (24)), because u is only added to the environment (sites 
outside of the impurity) while /i affects every site in the lattice, 
including the impurity. This difference allows us to use the 
two-step self-consistency scheme to determine /i, as shown 
in Fig. 1. Specifically, we first fit ji at the mean-field stage, 
to ensure (n) is correct. Then at the embedding stage, we 
vary ji and u simultaneously, following Eq. (32). This means 
that the DMET mean-field solution (and thus definition of the 
impurity model) stay the same, but the relative energy levels 
of the impurity change as compared to the bath states, which 
allows us to adjust the filling on the impurity. 

Fitting at the embedding stage means we need to solve the 
correlated impurity problem more than once in a single DMET 
self-consistency iteration. This increases the computational 
cost. Our strategy is to allow only one iteration of chemical 
potential fitting in each DMET iteration, corresponding to at 
most three DMRG calculations. Because fitting /i is a one 
dimensional search, even with this crude approach, we can 
usually control the relative deviation of (n) to under 10“^. 


II. ERROR MODEL 

As described in the main text, we consider 3 sources of 
error: (i) errors in DMET self-consistency, (ii) finite M in 
the DMRG solver, and (iii) finite impurity cluster size. The 
DMET self-consistency error is estimated as — 

where and are the energies of the last two 

DMET self-consistency iterations. A typical DMET calcula¬ 
tion oscillates between two slightly different solutions with 
the magnitude of the oscillations decreasing with the num¬ 
ber of iterations. We use the range of oscillation as a rep¬ 
resentation of the self-consistency error. The error distribu¬ 
tions across the range of calculations in this work are shown 
in Fig. 2, with the average values in the inset. For most 
points in the phase diagram, and for all cluster sizes, the 
self-consistency error is less than O.OOObt. For 4x4 clus¬ 
ters DMET calculations are the harder to converge, giving a 
largest error of up to 0.002t, and an average self-consistency 
error approximately twice as large as that for the other cluster 
shapes. 

For impurity clusters larger than the 2x2 cluster (where 
our DMRG solver essentially does an exact diagonalization), 
there is error due to using finite M in the DMRG impurity 
solver. The error due to finite M has two components: 

1. variational error in the DMRG calculation, which is 
usually assumed proportional to the density matrix trun¬ 
cation weight 

2. the DMET correlation potential error 6u, as Su is a func¬ 
tion of the impurity density matrices, and these have an 
error for finite M. 

For the 4 X 2 and 8x2 clusters, Su appears negligible. For 
these clusters, we carry out the DMET self-consistency with 
lower M to obtain the DMET correlation potential u, then do 
a few final DMRG calculations at large M to extrapolate to the 


180 
150 
120 
90 
60 
30 

0.0000 0.0005 0.0010 0.0015 0.0020 

FIG. 2. Distribution and average value (inset) of the DMET self- 
consistency error in the energy (units of t) for each cluster size. 



M ^ oc exact solver limit. For 4x4 clusters, the U = 2 data 
is processed in this way as well. However, for other values 
of U using the 4x4 clusters, the DMRG truncation weight 
is as large as 10“^ for low to intermediate doping with our 
accessible M, thus making the contribution of Su also signif¬ 
icant. To compensate for this, we first carry out the DMET 
self-consistency with a series of different M up to 1200, and 
linearly extrapolate the energy to the M = oc limit, Ei . This 
thus extrapolates errors from both source 1 and 2, assuming 
Su ^ Suj. Another further set of DMRG calculations are then 
done with M up to 2000, using the converged correlation po¬ 
tential from the DMET self-consistency with the largest M. 
This second set of results are then extrapolated again against 
the truncated weight to obtain an energy E 2 , which only ac¬ 
counts for the error from source 1. Although the linear relation 
between the source 2 error and the truncation weight need not 
hold in general, in practice, we find that Su = \\Ei — E 2 \ 
gives a crude estimate of Su^ Therefore, we report the 4 x 4 
cluster energy as E^ 4 x 4 = E E 2 ), with a final uncer¬ 
tainty of ’ where SEi is a combina¬ 

tion of the linear regression uncertainty and the uncertainties 
of the original data points (from DMET self-consistency er¬ 
ror), while E 2 does not have any self-consistency error. Fig. 3 
illustrates the set of computations and linear extrapolations 
performed with each 4x4 cluster to obtain the 4 x 4 cluster 
energy and error estimate. 

After obtaining the energy and observables for each clus¬ 
ter size, we extrapolate to the thermodynamic limit using the 
relation AE^^ oc Since both the 4 x 4 and 2x8 

clusters are 16 site clusters, we must choose which one to use 
in the extrapolation to the thermodynamic limit. We believe 
that 4x4 clusters have less finite size error than the 8 x 2 
clusters, and thus we generally use these in the extrapolation. 
However, at certain points in the phase diagram (e.g. at strong 
coupling, or negative t') there is a strong tendency to inhomo¬ 
geneity, and the 4 X 4 clusters cannot accommodate the new 
order parameters, resulting in a much higher energy. In such 
cases, namely, when (a) 4 x 4 and 8x2 clusters show different 
orders, and (b) the 8x2 cluster is lower in energy, we use the 
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(a)t/=4E =-1.033(2) 
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FIG. 4. Cluster size extrapolation fox U = 4, t' = 0 at various 
fillings. The black dots are finite size results. The red error bars are 
the confidence intervals for the thermodynamic limit. 


(b)C=6 E = -0.866(2) 



5 

(c)t/=8 E = -0.748(4) 


FIG. 3. Computations involved in the estimate of the 4x4 clus¬ 
ter DMET energy. The black dots (with error bars for the self- 
consistency error) are DMET self-consistent results using different 
DMRG M. These points are extrapolated to obtain Ei . The red dots 
are DMRG results using the “best” self-consistent correlation poten¬ 
tial, which are then extrapolated to obtain E 2 . The final 4x4 cluster 
DMET energies are reported as = |(^i + ^ 2 ). The plots are 
shown for = 0, n = 0.875 and (a) U=4 (b) U=6 (c) U=S. 



(a)2 X 2 


(b)4 X 2 


(c)4 X 4 


0.10 0.12 0.13 0.15 0.15 0.13 0.12 0.10 



0.16 0.12 0.08 0.04 0.03 0.07 0.12 0.16 

(d)8 X 2 


EIG. 5. Local order parameters for 6/ = 4,n = 0.875, t' = 0. The 
legend is the same as for Eig. 5 in the main text. 


8x2 cluster energy for the extrapolation. 

The cluster size extrapolation works surprisingly well given 
the limited number and small sizes of the clusters, although it 
contributes the main source of error in the final uncertainty. In 
Fig. 4 we show some of the extrapolation results at 7/ = 4. 
At half-filling and in the overdoped region (n < 0.8), the 
linear relation used in the cluster size extrapolation appears 
quite good even for these small clusters. In the underdoped 
region, however, the energy is more strongly dependent on 
the cluster shape, often because the system has a strong ten¬ 
dency to establish an inhomogenous phase. In Fig. 5, we plot 
the local order parameters at n = 0.875, where the 8 x 2 clus¬ 


ter calculation gives an incommensurate antiferromagnetic or¬ 
der. Although the 8 X 2 cluster energy (—1.0288) is slightly 
higher than the 4 x 4 cluster result (—1.033), its inhomogene¬ 
ity suggests the existence of a low-lying inhomogeneous state 
that can be (relatively) stabilized by special cluster shapes. 
Nonetheless, even in the underdoped region, the error model 
appears to give a reliable estimate of the energy at the thermo¬ 
dynamic limit, albeit with a large uncertainty. 

Fig. 6 shows the final energy errors for t' = ±0.2 across 
the phase diagram. The same plot for t' = 0 is shown in 
Fig. 2 in the main text. The overall uncertainty for = 0.2 is 
smaller than t' = 0 (see Fig. 2 in the main text) and t' = —0.2, 
as is the maximum uncertainty (0.011 compared to 0.037 and 
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Z) 


Z) 



(a)t' = 0.2 




FIG. 8. Staggered magnetization (m) of the half-filled Hubbard 
model for t' = ±0.2 and t' = 0. 



(b)F = -0.2 

FIG. 6. DMET energy uncertainty plot for the frustrated Hubbard 
model with t' = ±0.2. Refer to Fig. 2 in the main text for the 
legend. 



FIG. 7. Examples of thermodynamic extrapolations where the energy 
is sensitive to cluster shape. 

0.02t, respectively). As mentioned before, the main source of 
error is the cluster size extrapolation. Two examples of large 
uncertainties due to cluster size (and shape) effect are shown 
in Fig. 7. The largest uncertainties are observed at f/ = 6 and 
moderate doping. 

III. FURTHER RESULTS 

In this section, we will expand on the determination of the 
phase diagram (Fig. 3 in the main text). 

The staggered magnetization for the frustrated Hubbard 


FIG. 9. At U = 2, the uncertainty of antiferromagnetic order param¬ 
eter decreases exponentially with doping. The exponent is 65 ± 4. 


model at half-filling (compared to the t' = 0 model and the 
Heisenberg limit) is shown in Fig. 8. Due to particle-hole 
symmetry, the plot is identical for t' = ±0.2. The onset 
of antiferromagnetism is at finite U in the frustrated model, 
between U = 2 and 3.5, consistent with previous quantum 
Monte Carlo simulations^^. The large error bar at f/ = 3 indi¬ 
cates the sensitivity to impurity cluster sizes near the phase 
boundary, resulting in a large uncertainty in the thermody¬ 
namic extrapolation. 

At weak coupling U = 2, we find that the antiferromag¬ 
netism (in the non-frustrated model) is destroyed already at 
small doping x = 0.05, where the staggered magnetization is 
m = 0.00 ± 0.05. Although the expectation value is 0, the 
relatively large uncertainty 6m reflects that short-range spin 
fluctations are still signiflcant, although long-range order does 
not exist. As we increase doping, Sm decreases exponentially 
(Fig. 9). At (7 = 2, we do not And d-wave superconductivity, 
to within numerical precision. 


0.11 0.11 0.13 0.15 0.15 0.13 0.11 0.11 



FIG. 10. Inhomogeneous order from 8x2 cluster calculations at 
U = 4,t' = -0.2 n = 0.875. 
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TABLE L Energy comparison for different 16-site impurity clusters 
at f/ = 4 and t' — —0.2. 


n 

Esx2 

EaxA 

0.8 

-1.10483(6) 

-1.0507(4) 

0.85 

-1.0162(1) 

-1.020(2) 

0.875 

-0.9966(1) 

-0.9989(7) 


We now discuss U = 4. We have already shown the 
order parameters, and the observed thermodynamic extrapo¬ 
lated ground state orders are all homogeneous. However, for 
t' = —0.2, the 8 X 2 cluster calculations result in an inhomo¬ 
geneous state at doping n = 0.8 — 0.875, although the energy 
is significantly higher than obtained with the 4x4 clusters at 
the same fillings. An example of an inhomogeneous pattern is 
shown in Fig. 10, where one can see a pair density wave and 
incommensurate magnetic order. In Table I, we compare the 
energies between the 8 x 2 cluster and 4x4 cluster results at 
relevant points in the phase diagram for f/ = 4. In all these 
cases, the 8 X 2 cluster has a higher energy, suggesting that 
the ground state at f/ = 4 is homogeneous or inhomogeneous 
with a very long wavelength that does not fit in our cluster 
shapes. 


0.04 0.08 0.14 0.19 0.23 0.15 0.11 0.05 



0.34 0.28 0.19 0.03 0.16 0.16 0.23 0.32 


(a)n = 0.875 

0.07 0.12 0.17 0.20 0.22 0.17 0.14 0.08 


0.30 0.22 0.11 0.06 0.16 0.09 0.18 0.29 

(b)n = 0.85 


0.21 0.20 0.19 0.19 0.19 0.20 0.20 0.21 



0.02 0.00 0.04 0.06 0.06 0.03 0.00 0.02 

{c)n = 0.8 


EIG. 11. Inhomogeneous order from 8x2 cluster calculations at 
U — Q and = 0 with fillings 0.875 to 0.8. 

AiU = 6, more interesting inhomogeneous orders start 
to appear. At = 0, 8 x 2 clusters result in various orders 
(Fig. 11). At both n = 0.875 and n = 0.85, 4x4 clusters are 







llall 

IMi 

liJ| 

, t 



t=f 


FIG. 12. Local order parameters from 4x4 cluster calculations at 
U — Q and = 0 with fillings 0.875 and 0.85. 


significantly lower in energy, suggesting the charge, spin and 
pairing orders shown in Fig. 11(a) and 11(b) are not stable. 
Atn = 0.875, a homogeneous solution with both supercon¬ 
ductivity and antiferromagnetism is found (Fig. 12(a)). How¬ 
ever, the thermodynamic extrapolation gives zero for both AF 
and SC order parameters. At n = 0.85, the 4 x 4 cluster 
result also shows slight inhomogeneity, with a (tt, tt) modula¬ 
tion of the d-wave order parameter (Fig. 12(b)). At n = 0.8, 
where the 8x2 impurity cluster gives a slightly lower energy 
{AE = 0.003(2)), where DMET calculations indicate a weak 
spin density wave (Fig. 11(c)). This spin density wave may 
still exist in the thermodynamic limit because the amplitude is 
comparable to the staggered magnetization in smaller clusters 
(eg. m = 0.04 for 2 x 2 clusters). 



We now turn to t' = —0.2. At n = 0.8 and 0.875, 
8x2 cluster calculations show inhomogeneous orders. At 
n = 0.875, the pattern is similar to what we observed for 
t' = 0 at the same filling, and its energy Esx 2 = —0.8402(4) 
is much higher than that of the 4 x 4 homogeneous solution 
F^ 4 x 4 = —0.850(3). At n = 0.8 (Fig. 13), both 4x4 and 
8x2 cluster calculations show tt -phase shifts in the spin den¬ 
sity and d-wave order, while the 8 x 2 cluster has an additional 
charge density wave. They are very similar in energy, with 
^ 8 x 2 = -0.9283(2) andF^ 4 x 4 = -0.927(3). This suggests 
that the ground state here is superconducting with a superim¬ 
posed spin density wave. 

Most results for the underdoped region at f/ = 8 are al¬ 
ready shown in the main text (Fig. 5). In Table II, we compare 
energies for the two 16-site clusters. At all the points shown 
in the table, the 8 x 2 cluster gives a lower energy. An unusual 
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result is that at n = 0.8, = 0, the 8x2 cluster shows a 

homogeneous solution, while both the 4x4 and 4x2 clusters 
give a spin 7r-phase shift. This unusual behaviour, where the 
8x2 solution favours homogeneity, is related to the large error 
in the energy at this point in the thermodynamic extrapolation 
{dE = 0.03). 


TABLE II. Energy comparison for different 16-site impurity clusters 
at f/ = 8. 


t' 

n 

E8x2 

±4X4 

0 

0.8 

-0.9018(13) 

-0.873(6) 

0 

0.875 

-0.7548(4) 

-0.748(4) 

-0.2 

0.8 

-0.8487(4) 

-0.846(10) 

-0.2 

0.875 

-0.7556(5) 

-0.737(7)" 


^ The error estimate may not be reliable at this point, because we have only 
two self-consistent DMET calculations with M=1000 and 1200. 


0.19 0.22 0.21 0.19 0.20 0.22 0.20 0.17 



0.17 0.04 0.09 0.17 0.14 0.03 0.10 0.20 

(a)L = 4, AE; = 0.0024 ± 0.0005 


0.17 0.23 0.23 0.17 0.17 0.23 0.23 0.16 



0.26 0.09 0.10 0.24 0.23 0.09 0.10 0.26 

(b)U = 6,AE = -0.001 ± 0.003 


0.15 0.25 0.24 0.16 0.20 0.26 0.22 0.13 



0.32 0.09 0.15 0.29 0.22 0.03 0.17 0.34 

(c)U = 8,AE = -0.03 ± 0.010 


EIG. 14. Local order parameters forn = 0.8, = —0.2 with differ¬ 

ent U/t. AE = E 8 x 2 — L^ 4 x 4 is the energy difference between the 
two 16-site clusters. 

In the above we discussed competing orders at different 
coupling strength, but it is also interesting to look at their evo¬ 
lution with U, as shown in Fig. 14. When U increases from 
4 to 8, we see increasing charge and spin inhomogeneity al¬ 
though they all show the same pattern of charge localization 
and a spin tt- phase shift. The d-wave pairing strength first 
increases and then becomes inhomogeneous. The energy dif¬ 
ference between the 8 x 2 and 4x4 clusters also changes 


monotonically, although the large error bars prevent us from 
definitively determining the true order at f/ =6 and 8. 

Finally, we end our discussion on the results by showing 
the energies across the phase space in Fig. 15. At half-filling, 
the energy in the frustrated model t' = ±0.2 is slightly below 

= 0, while the difference becomes negligible at large U. At 
the large doping, eg. n < 0.8, the energy order is dominated 
by the kinetic effects, i.e. F^t'=-o .2 > > ^t'=o. 2 - The 

energy curves show more complicate characteristics at under¬ 
doped region, especially for = 0 and t' = —0.2. 



(a) 



(b) 



(c) 



n 


(d) 

FIG. 15. DMET thermodynamic energy over the phase space, (a) 
U=2 (b) U=4 (c) U=6 (d) U=8 
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IV. DATA SET 



FIG. 16. The encoding of local order parameters for all impurity 
clusters. Numbers shown in the circles represent the order of sites, 
which is associated with charge density and spin density. The num¬ 
bers in the rhombus represent the order of bonds, or pairs between 
neighbor sites, which is associated pairing strength. Some numbers 
are omitted since it is easy to know what they are. 


In the attached TDL.csv file, we present the energy, chem¬ 
ical potential and (averaged) order parameters computed and 
their uncertainties at the thermodynamic limit. Since the aver¬ 
aged order parameters are meaningless when inhomogeneity 
dominates, we have removed these entries from the table. 


In the file clusters.csv, we present the results for finite im¬ 
purity clusters. In addition to the results available at thermo¬ 
dynamic limit, we also present the local order parameters. The 
local order parameters are encoded in an ID array, which is 
explained in Fig. The errors shown only include the DMET 
convergence error, as the other sources of error can be de¬ 
duced using the procedures described above, from the raw 
data. We also include the local orders (charge, spin and pair¬ 
ing strength) in this table as a ID array. The order of the sites 
and pairs are shown in Fig. 16. 

V. APPENDIX 

Here we prove the equivalence of the Fock spaces spanned 
by Cl and C 2 in the construction of the impurity Schmidt sub¬ 
space, as defined in section (IC). Precisely, we need to prove 

1. C 2 is orthonormal, CJC 2 = I. (It is easy to see Ci 
is orthonormal, because Q is a unitary matrix from QR 
decomposition). 

2. C 2 = CiV, which is equivalent to C 1 C 2 = V, where 
V is unitary. 

To prove (1) C 2 C 2 = L we need the idempotency of den¬ 
sity matrices = p. Considering only the upper-left block 
of p, we have 

Amp Pc Pc — Pimp (33) 

From Eq. (13) and (14), we know A' = U{I — A)^, B' = 
—pc{A'^)~^. Therefore 


( A^ A-^pI \ f A A' \ _ f A^A^A-^pIp^{A^)-^ A^A' - A-^pip^{A'^)-^ 

A'^A' + {A')-^pIpM'^)-^ 

f A +A-U(/-A)A-5 A5(/-A)5 - A-U(/-A)(/-A)-5 \ _ 

V (/-A)5A5 - (/-A)-5A(/-A)A-5 /-A + (/-A)-5A(/-A)(/-A)-5 / 

(34) 


For (2), since 


V = CfC2 = 


A 

Q'^Pc{A'^)~^ 


A' \ 

-Q'^ Pc{A''^) 1 ) 


we have 


A 

A' \ 

( A^ A-^p^Q \ 

Q'^Pc{A'^)~^ 

-Q^Pc(A'^) ^ ) 

1 

T 


pTQ _ pTQ 

Pc - Pc Pc(.AA^)~^ p'^Q + Q'^pc{A'A''^)~^p'^Q 


I 0 \ 

0 + 


(35) 


(36) 


In the bottom-right block 

A-i + (/ - A)-i =UA-^U'^ + U{I - A)-if/^ 
=f/A-i(/-A)-if/^ 

=[A(/-A)]-i 

= (R^R)-^ = R-\R^)-^ 
(37) 


So VV^ = I. Here we assume R is invertible, which is true 
if and only if we have the full set of Nc bath orbitals coupled 
to the impurity. This is generally true in lattice settings where 
the impurity and the environment are strongly coupled. Some- 
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times the bath can be smaller than the impurity in molecules and when we use a large basis set, and in these cases, special 

treatment is needed. 
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